function [ sensors ] = espass_alg( dataset, k, m, e )
%espass_alg eSPASS algorith described in section 4, alg 4
%   Intelligently picks placement and schedule, optimizes for worst case.
%   Picks m sensors from the dataset and distributes across k times.
%   e is the threshold after which the algorithm is happy with the solution

    % Initialize values
    cmin = 0;
    v = 1;
    for i=2:size(dataset, 2)
        v = union(v, i);
    end
    cmax = get_quality(dataset, v, Inf);
    beta = 68/96;
    abest = [];
    
    % Loop until cmax and cmin are close enough to stop, based on threshold
    while((cmax - cmin) >= e)

        c = (cmax + cmin) / 2;
        
        % Compute the big set, by looping over all sensors
        big = [];
        for s = 1:size(dataset, 2)
            if(get_quality(dataset, s, c) >= (beta * c))
                big = union(big, s);
            end
        end
        kprime = k;
        
        % Ats matrix stores all of the elements and their time slots
        ats = zeros(k, m);
        % Index of first zero in each time slot in the above matrix
        ats_zero = ones(k, 1);
        % Loop over each element of big
        for s = 1:size(big, 2)
            % Add s to kprime slot
            ats(kprime, ats_zero(kprime)) = big(s);
            ats_zero(kprime) = ats_zero(kprime) + 1;
            kprime = kprime - 1;
            
            % If kprime = 0, then we need to get more optimistic
            if (kprime == 0)
                cmin = c;
                abest = ats;
                abest_zero = ats_zero;
                break;
            end
        end
        % Move on to next while iteration, if kprime == 0
        if (kprime == 0)
            continue;
        end
        
        % Initialize prime values for next loop
        vprime = setdiff(v, big);
        mprime = m - size(big, 2);
        
        % Run GAPS for time slots up to kprime
        % First have to create a dataset with only vprime columns
        dataprime = [];
        for i=1:size(dataset, 2)
            if (ismember(i, vprime))
                dataprime = [dataprime, dataset(:, i)];
            end
        end
        ssmall = gaps_alg(dataprime, kprime, mprime, c);
        % Convert the gaps output into standard ats
        asmall = zeros(kprime, m);
        asmall_zero = ones(kprime, 1);
        for s = 1:size(ssmall, 2)
            asmall(ssmall(2, s), asmall_zero(ssmall(2, s))) = ssmall(1, s);
            asmall_zero(ssmall(2, s)) = asmall_zero(ssmall(2, s)) + 1;
        end

        % Calculate the sum of all time slots
        qualsum = 0;
        for t = 1:kprime
            % From here on use dataprime, instead of dataset. Or should we?
            qualsum = qualsum + get_quality(dataprime, asmall(t, :), Inf);
        end
        
        % If this value is small, we overdid things, so become pesimisitic
        if(qualsum < (kprime * c / 2))
            cmax = c;
            continue;
        end
        
        % Loop over all sensor combinations <= kprime
        for i = 1:kprime
            for j = 1:kprime
                % Check all the continue condtions
                if (i == j)
                    continue;
                end
                % Check the other conditions, and end loop iteration if met
                if (get_quality(dataprime, asmall(j, :), c) > (beta * c))
                    continue;
                end
                % Probably don't have to do this one for every j
                if (get_quality(dataprime, asmall(i, :), c) < (3 * beta * c))
                    continue;
                end
                
                % Loop over all elements in i, and move them to j until j
                % becomes better than beta * c
                % This works because of silly static matlab for loops
                for s = (asmall_zero(i) - 1):1
                    % Move the element from i to j
                    asmall(j, asmall_zero(j)) = asmall(i, s);
                    asmall_zero(j) = asmall_zero(j) + 1;
                    asmall(i, s) = 0;
                    asmall_zero(i) = asmall_zero(i) - 1;
                    
                    % Check to see if j is quality enough to stop
                    if (get_quality(dataprime, asmall(j, :), c) >= (beta * c))
                        break;
                    end
                end
            end
        end
        
        % If we got to this point, then get more optimistic
        cmin = c;
        % Have to copy elements from ats above and asmall into abest
        abest = ats;
        abest_zero = ats_zero;

        % Copy and map entries in asmall back to sensor IDs from V
        for t = 1:kprime
            for i = 1:(asmall_zero(t) - 1)
                abest(t, abest_zero(t)) = vprime(1, asmall(t, i));
                abest_zero(t) = abest_zero(t) + 1;
            end
        end
    end
    
    % Build the list of sensors from abest
    sensors = zeros(2, m);
    nextsensor = 1;
    % Loop over each row of ats
    for t = 1:k
        for i = 1:(abest_zero(t) - 1)
            sensors(1, nextsensor) = abest(t, i);
            sensors(2, nextsensor) = t;
            nextsensor = nextsensor + 1;
        end
    end

end

